The Cross- Wavelet Transform and Analysis of Quasiperiodic 
Behavior in the Pearson-Readhead VLBI Survey Sources 
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ABSTRACT 

We introduce an algorithm for applying a cross- wavelet transform to analysis of quasiperiodic 
variations in a time-series, and introduce significance tests for the technique. We apply a contin- 
uous wavelet transform and the cross-wavelet algorithm to the Pearson-Readhead VLBI survey 
sources using data obtained from the University of Michigan 26-m parabloid at observing fre- 
quencies of 14.5, 8.0, and 4.8 GHz. Thirty of the sixty-two sources were chosen to have sufficient 
data for analysis, having at least 100 data points for a given time-series. Of these thirty sources, 
a little more than half exhibited evidence for quasiperiodic behavior in at least one observing 
frequency, with a mean characteristic period of 2.4 yr and standard deviation of 1.3 yr. We find 
that out of the thirty sources, there were about four time scales for every ten time series, and 
about half of those sources showing quasiperiodic behavior repeated the behavior in at least one 
other observing frequency. 

Subject headings: galaxies: active — methods: data analysis — radio continuum: galaxies 



1. INTRODUCTION 

It is well accepted that centimeter waveband 
emission from AGNs is associated with a jet of 
synchotron plasma, the accretion structure and 
immediate environment of a supermassive black 
hole contributing broad band emission from the 
infrared to the gamma ray spectrum. Temporal 
variations are observed in the radio flux, and a 
number of processes have been proposed to explain 
this, such as an accretion rate that may change 
with time, an accretion disk that exhibits insta- 
bility, an outflow that may be Kelvin-Helmholtz 
unstable, or an outflow that may interact with am- 
bient inhomogeneities. The temporal variations in 
the radio flux motivate a search for characteristic 
time scales, and, if found, would lend insight into 
the mechanisms causing the variations. 

In practice, searching for possible quasiperiodic 
behavior masked by a stochastic component can 
be rather difficult, especially in the context of ir- 
regular time sampling. One promising technique 
is to perform a continuous wavelet transform on 



the signal, and map out the coefficients in wavelet 
space. This was done on data for the BL Lac OJ 
287, resulting in evidence for periods in the radio 
spectral region of ~ 1.66 yr, and ~ 1.12 yr domi- 
nating in the 1980s (Hughes, Aller, & Aller 1998). 
The results led Hughes, Aller, & Aller (1998) to 
propose that the results can best be explained by 
a "shock- in- jet" model, in which the longer peri- 
odicity is linked with an otherwise quiescent jet, 
and the shorter with a passing shock-wave. 

The continuous wavelet transform detects 
quasiperiodic behavior through visual examina- 
tion of the map over translation and dilation, and 
assigns a characteristic time scale by looking for 
peaks in the time-averaged wavelet power (see 
§ 4.1). The scales that the peaks occur at cor- 
respond to Fourier periods, from which a charac- 
teristic time scale may be deduced. Although the 
continuous transform is effective, when assigning 
a characteristic time scale it does not fully include 
information of how a time series varies in dilation; 
the time scale corresponds to those dilations where 
the time-averaged wavelet power spectrum peaks. 



1 



In order to use all of the available information in 
the transform, we examine the cross- wavelet trans- 
form, in which the coefficients for the continuous 
transform of a signal arc multiplied by the complex 
conjugate of the coefficients of another signal. The 
results are then mapped out in wavelet space and 
analyzed for correlation. A cross-wavelet analysis 
was used on the OJ-94 project light curve (Lehto 
1999), and has seen application in other areas of 
science as well (Pancheva et al. 2000; Kyprianou 
& Staszewski 1999). We examine this method in 
greater depth, as well as a technique for finding 
the characteristic period of a signal and apply 
this to the Pearson-Readhead (PR) VLBI survey 
sources using data from the University of Michi- 
gan Radio Astronomy Observatory (UMRAO) ob- 
served at three frequencies. The PR group is well 
suited for analyzing quasiperiodicity, due to the 
high signal-to-noise in the data for most of these 
sources (Aller, Aller, & Hughes 2002), as well as 
investigating quasiperiodic variations (QPVs) as 
a function of optical source type, since it contains 
representatives from optical classes QSO, galaxy, 
and BL Lac. In addition, the time base of the 
UMRAO data is around twenty years for most 
sources, allowing possible quasiperiodic events of 
as long as four years to be observed for at least 
four to five cycles. Using UMRAO data, evidence 
has been found for periodic behavior in the cen- 
timeter waveband emission for the BL Lac object 
0235+164 (Roy et al. 2000). However, there have 
not been any convincing indications of periodicity 
for PR sources (Aller, Aller, & Hughes 2002). 

2. THE DATA 

The data for this study were acquired as part 
of the UMRAO program, using the University of 
Michigan 26-meter parabloid. Sixty-two of the 
sixty-five PR sources are within the declination 
limits of the parabloid, and these sources have 
been observed at least trimonthly at 14.5, 8.0, and 
4.8 GHz since the fall of 1984, with occasional gaps 
due to poor weather and positions too near to the 
sun to be observed. The observing technique and 
reduction procedures used are those described in 
Aller et al. (1985) and Aller et al. (2002), with the 
latter including a rescaling of the data for the pur- 
pose of conforming with the flux system of Ott et 
al. (1994). 



3. WAVELET TRANSFORMS 

3.1. The Continuous Wavelet Transform 

The continuous wavelet transform involves de- 
composing a signal f(t) into a number of trans- 
lated and dilated wavelets. The main idea behind 
this is to take a "mother" wavelet, translate and 
dilate it, convolve it with the function of interest, 
and map out the coefficients in "wavelet space" , 
spanned by translation and dilation. Periodic be- 
havior then shows up as a pattern spanning all 
translations at a given dilation, and this redun- 
dancy in the wavelet space makes detection of peri- 
odic behavior rather easy. The wavelet transform 
preserves temporal locality, which is an advantage 
over Fourier analysis. For instance, power associ- 
ated with irregular sampling does not contribute 
to the coefficients as in Fourier analysis, which 
is extremely helpful when using poorly-sampled 
data. 

There arc several common types of mother 
wavelets. In this analysis we use a Morlet wavelet 
of k^f, = 6, given by 

with coefficients 

/(M') = / f(t)^ t ,(t)dt (2) 

JR 

and 



V;t' (t) = 




ieR + , ten. (3) 



The continuous transform must also satisfy an ad- 
missibility condition, requiring zero mean: 

f iP(t)dt = 0. (4) 

JR 

It should be noted that I corresponds to dilations, 
and t' refers to translations. Generally the pa- 
rameter k^, remains fixed throughout the analysis, 
and for simplicity will hereafter be just k. Fig- 
ure la shows the continuous wavelet transform for 
a sinusoidal signal. The periodic behavior is eas- 
ily revealed by the redundancy in the plot. For 
the interested reader, more detailed information 
on the continuous wavelet transform may be found 
in Farge (1992). 
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3.2. The Cross- Wavelet Transform 



Equation (9) then becomes 



The continuous wavelet transform is effective 
for examining how a time series varies in time and 
scale, but fails to include how it varies over a range 
of scales when assigning a period that best charac- 
terizes it. After identifying that a periodic pattern 
exists in what can be potentially noisy and poorly 
sampled data, one finds the dilation which char- 
acterizes the period from the time-averaged data 
(i.e., the wavelet power spectrum), and from that 
dilation the period is calculated. For a quasiperi- 
odic signal there is no unique dilation; there is a 
need for a method of directly measuring a char- 
acteristic time scale, or scales, that includes infor- 
mation on how a source varies in dilation, and for 
this we shall examine the cross- wavelet transform. 

Given two time series f a (t) and f m (t), we can 
construct the cross- wavelet transform 

f c (l,t') = f a (l,t')-f m (l,t'), (5) 

where f a and f m are given by Equation (2). In 
this analysis f a {t) is the source signal and f m {t) 
is a sinusoidal mock signal, although one can use 
any two time series thought to be correlated. 

Because we are looking for quasiperiodic behav- 
ior, first assume an ideal signal of the form 

f a {t)=A a e i ^ t+ ^\ -oo<i<oo, (6) 

where to a = 27r/r a for period r a . Anticipating that 
it will be productive to cross the actual signal with 
a signal of a similar form, we try a mock signal 

f m (t) = A m e^ M ™\ -^<t<^. (7) 

Each of these signals will have continuous wavelet 
coefficients given by Equation (2): 

f(l,t') = W2Z 7 r 1 /2 e -^-fc)yM'+«. (8) 

The cross- wavelet (Eq. [2] ) for these two signals be- 
comes 



f c (l,t') = Ale-^-^e^'+^l 



(10) 



f c (l,t')=2^A a A r , 



-^[(uil+uj'^l 2 -2kl(oj a +u> m )] Mo a -oj m )it' 



(9) 



Defining 
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2^A a A m e 

W a + LO m 



— ^ a 
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which is easily interpreted. From this equation 
we note two important properties of this cross- 
wavelet — first, it has the form of a Gaussian in 
the dilation coordinate I, and second, it is sinu- 
soidal in the translation coordinate t' with fre- 
quency given by the difference in the frequencies 
of the actual and the mock signal (i.e., the beat 
frequency). When these two frequencies are equal, 
the translation dependence is lost, and the cross- 
wavelet reduces to a Gaussian in the dilation co- 
ordinate. Figure lb, shows the cross- wavelet for 
two sinusoids oscillating with the nearly same fre- 
quency 

3.3. Using the Cross- Wavelet Transform 
to Analyze Periodicity 

Motivated by the results from the previous sec- 
tion, we can develop an algorithm to find the pe- 
riod which best characterizes a signal, i.e., the 
characteristic time scale. Because we are dealing 
with real signals, noise will be present, and we ex- 
plore the role of noise in § 4 and § 5. In addition, 
the values of t and t' do not extend to infinity, 
but over the range t = a to t = b. Because of the 
finiteness of the range of t, the range of dilations 
be iconics I'min *C ^ <C Irnax • ^Ve can still use Equa- 
tion (8) as an approximation when our range of I 
and t' is such that the wavelets t^w (t) contribute 
negligibly outside of a < t < b, and 

fb r poo r 

/ e -[( u -f)«±*|~-f ](ft / e _ [(w _*) it± §|tM^ 

J a J — oo 

(ii) 

Returning to Equation (10), we define the mod- 
ulus as 



f c (l,t') =A 2 l 2 e-W 2 - 2k "< l \ (12) 



Integrating over the region of wavelet space we 
mapped out, weighting by the scale, and approxi- 
mating the integral over I out to infinity, we arrive 
at 



F c (r m ) 



fc(l,t') 



A 2 (b-a) 



I 



-dt'dl 



v V v 



(13) 
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We have empirically found that this function at- 
tains it maximum value at r a = 0.973r m for a 
value of k = 6. This result is not surprising, as 
one would expect the power of the cross coeffi- 
cients to reach their maximum when the periods 
for the two signals are almost equal. One can think 
of the mock signal as acting like a filter in wavelet 
space; however, the power of the mock coefficients, 
|/ m | 2 , are slightly asymmetric about the peak (see 
Eq. [8] ) . The weighted total power of the cross co- 
efficients, F c , does not attain its maximum exactly 
at T m = r a because of this slight asymmetry. 

To avoid confusion, one must remember that 
previously when the integration was over [a, b] we 
integrated over the original coordinate t and, be- 
cause the wavelets have compact support in t, we 
approximated the integral by taking the limits to 
be [—00,00]. Here, however, the integral is over 
the translation coordinate t', so we must restrict 
the limits to [a, b]. In addition, the approximation 
of the integral over I is valid so long as f c (l,t')/l 
does not peak in the dilation coordinate too close 
to I 

min Or Ijnax- For OUT purpOSCS, Imin 

= 25t and 

Imax = N data St/3, where N data is the number of 
data points and St is the time spacing. Typical val- 
ues of b — a for the PR Survey Sources are roughly 
20 years, and we require at least Ndata — 100 data 
points for a reliable cross-wavelet analysis, which 
results in < l min < 0.4 yr and l max < 12 yr. 
We search for periods within 0.5 < r < (b — a) /A 
years. We have inspected |/ c | 2 /' graphically to 
assure that the approximation is valid for periods 
T a > 0.7 yr for all sources having more than 100 
data points. 

In practice, one must worry about edge effects 
as the time series is finite. Because the convo- 
lution of Equation (2) introduces power into the 
wavelet coefficients / from the discontinuity at the 
edges of the time series, there is a region where 
the wavelet coefficients will be contaminated from 
edge effects. This region is known as the cone of 
influence. The wavelet power associated with edge 
effects becomes negligible for translations t' far- 
ther than y/2l from the edge (Torrence & Compo 
1998) 

One more useful quantity is the point I where 
f c (l,t') peaks in the dilation coordinate. This is 



given by 

~ k K + u m ) + V(fc 2 + 4)(w2 + ujI) + 2k 2 U a Ur, 



When cj m = cj a , this becomes 



I 



fc + VP~+~2 

2uj a 



or, 



AttI 



2tt 



(14) 



(15) 



(16) 



k + ^/W+2 UJ a 

This equation shows the relationship between a 
scale I and the corresponding Fourier period, and 
can also be arrived at by inputing a sinusoid for 
f(t) in Equation (2). For clarification, we will use 
the term 'Fourier period' when referring to the 
period associated with a certain scale I, and the 
terms 'analyzing period' or 'mock period' when re- 
ferring to the period r m associated with the mock 
signal. 

4. SIGNIFICANCE TESTS 

4.1. Tests for the Continuous Transform 

Wavelet analysis has historically suffered from a 
lack of statistical significance tests. In Torrence & 
Compo (1998) an excellent discussion of statistical 
significance for the continuous wavelet transform 
is given and supported by Monte Carlo results, 
and we summarize the relevant points. When im- 
plementing the wavelet transform, one must use 
sums and discrete points rather than the theoreti- 
cal treatment given earlier, and we switch to using 
these. 

The model of correlated noise most likely to 
closely resemble the UMRAO data is the univari- 
ate lag-1 autoregressive (AR(1)) process (Hughes, 
Aller, & Aller 1992), given by 



-1 + z n , 



(17) 



where a is the assumed lag-1 autocorrelation and 
z n is a random deviate taken from white noise 
(note that this is a 'first order' process). The nor- 
malized discrete Fourier power spectrum of this 
model is 



I -a' 



1 + a 2 - 2acQs{2ir5t/T j y 



(18) 
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where Tj is the associated Fourier period (Eq.[16]) 
of a scale lj. If the time series of interest is an 
AR(1) process, then a slice in the modulus of the 
continuous coefficients, \f(l,t')\ 2 , along dilation 
at a constant translation (the local wavelet power 
spectrum) should have have the form in Equation 
(18). Note that a white noise process is given by 
a = 0. 

Using the theoretical background spectrum 
given above, we can develop significance tests for 
the continuous transform. The background spec- 
trum is the mean power spectrum expected for the 
assumed noise process (i.e., the 'background' pro- 
cess), against which we wish to compare the actual 
signal. Physical processes that are the result of 
this background process will produce power spec- 
tra that are normally distributed about this mean 
background spectrum. Assuming that the values 
in our time series f(t) are normally distributed, 
and because the square of a normally distributed 
variable is chi-square distributed with one degree 
of freedom, we expect the wavelet power |/| 2 to 
be chi-square distributed with two degrees of free- 
dom. The additional degree of freedom comes 
from the fact that the wavelet coefficients / are 
complex for the Morlet wavelet. In addition, the 
expectation value for the wavelet power of a white 
noise time series is just the variance a 2 ; this expec- 
tation value provides a convenient normalization. 
Using this normalization, the distribution for the 
local wavelet power spectrum is 



pXv_ 

=>±-J v , 



(19) 



been shown that the global wavelet spectrum pro- 
vides an efficient estimation of the true power of a 
time series (Percival 1995). Averaging the wavelet 
power spectrum as in Equation (20) increases the 
significance of the peaks, as the degrees of free- 
dom is increased beyond what is used in the local 
wavelet power spectrum. However, because the 
coefficients are correlated in both time and scale, 
the degrees of freedom for the global wavelet spec- 
trum are 



Vj =2W1 + 



NjSt 
ljSt 



(21) 



where Sto describes the decorrelation length in 
time. For the Morlet wavelet of k — 6, <5to = 2.32. 
The global wavelet spectrum of a sinusoid is shown 
in Figure 2. 

Alternatively, we can smooth the wavelet spec- 
trum in scale. We define the scale-averaged 
wavelet power as 



0=31 



(22) 



The scale-averaged wavelet power can be inter- 
preted as a time series of the average variance in 
the band lj 1 < lj < lj 2 . This distribution can be 
modeled as 



flit'i) 



(23) 



where the scale-averaged theoretical spectrum P 
is 



P= > ^. 



(24) 



where '=>' means 'is distributed as' and v is the 
degrees of freedom, in this case two. The indices 
on the scale I run from j = 1, 2, . . . , J, where J 
is the number of scales, and the indices on the 
translation t' run from 2 = 1,2,..., Ndata- 

We can also define the time-averaged wavelet 
power spectrum, or the global wavelet spectrum, 
as 



The degrees of freedom v\ in Equation (24) are 
modeled as 



where 



n = — 9 -\l 

J^mid V 



N t Sj 
Sjo ) 



(25) 



1 



(20) 



where ij and i'j are the indices of the initial and 
final translations t\ outside of the cone of influ- 
ence at a given scale lj, and Nj is the number t\ 
outside the cone of influence at that scale. It has 



^ h 



J, ., — 1 . <?0.5(j 1 +j 2 )6j 

I ■ — 1 oO'-!)^ 1 — 19 T 
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Ni is the number of dilations averaged over, and 
Sjo is the decorrelation distance in scale. For the 
Morlet wavelet of k = 6, 5j = 0.60. The 5j de- 
scribes how the dilations lj are constructed; it is 
common to construct them as fractional powers 
of two as it allows more emphasis on the smaller 
scales. The factor L avg /L m id corrects for the loss 
of degrees of freedom from dividing the wavelet 
power spectrum by scale in Equation (22). 

4.2. Tests for the Cross Transform 

Using the results of Torrence & Compo (1998) 
from the previous section, we can derive signifi- 
cance tests for the cross-wavelet transform. The 
power of the cross- wavelet f c can be written as 



(26) 



and is distributed as 



fa(lj,ti) fmiljttj) 



(7 



fm(lj, ij) 



2 Y 2 



(27) 

where the degrees of freedom v is two. We can also 
define the cross-wavelet global power spectrum, 



1 3 i 



(28) 



Noting that the power of the continuous coeffi- 
cients for the type of mock signal given in Equa- 
tion (7) is independent of the translation t', the 
time-averaged form of Equation (13), F c , then be- 
comes 



F c (r m ) = x: 



J P cg{1j) ^ P AG {h) fm{lj) 



J t2 



(29) 

where f\a(h) IS the global wavelet power spec- 
trum for the actual signal f a (ti)- This distribution 
can be modeled as 



for Vj and v\ given in the previous section. This 
allows us to smooth the wavelet power in both 
time and scale, thus increasing the degrees of free- 
dom, and to receive information of the charac- 
teristic time scale(s) of our time series from the 
mock coefficients f m rather than from the global 
wavelet spectrum. As mentioned in § 3.3, Equa- 
tion (29) allows us to interpret the cross-wavelet 
of a time series and an analyzing signal as a fil- 
ter in wavelet space, where the coefficients of the 
continuous transform for the analyzing signal fil- 
ter the global wavelet spectrum of our time series, 
f\ G - Averaging the cross transform in scale then 
allows us to find the analyzing signal with which 
the time series is best correlated over the entire 
range of time and scale. A plot of F c for a sinu- 
soid can be seen in Figure 2. 

' The results of this section are only valid when 
the cross- wavelet consists of one source that has an 
assumed background spectrum that is chi-square 
distributed. If one is crossing two sources, fi(t) 
and f2(t), that have assumed background spectra, 
Pi and P 2 , that are chi-square distributed, then 
the cross- wavelet is distributed as 



(71(72 



(32) 



where o\ and 02 are the respective standard de- 
viations, the degrees of freedom v is two for com- 
plex wavelets, and Z v (p) is the confidence level 
for a given probability p for the square root of the 
product of two chi-square distributions (Torrence 
& Compo 1998). One can find the confidence level 
Z v by inverting the integral p = J Q " f v {z)dz. The 
probability distribution is given by 



fM = 



->2-v 



I>/2) 



K (z), 



(33) 



where z is the random variable, T is the Gamma 
function, and K (z) is the modified Bessel func- 
tion of order zero. 



F c (r m ) 



J Pj \fm{h) 
3 = 1 



-f > (30) 



V 3 



where the degrees of freedom u'j are modeled as 
v'^Vj + vu (31) 
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5. ALGORITHM FOR IMPLEMENT- 
ING THE CROSS TRANSFORM TO 
ANALYZE QUASIPERIODIC BEHAV- 
IOR AND SIMULATIONS 

5.1. Algorithm for Analysis of Quasiperi- 
odic Behavior 

When using real data, we have no a priori 
knowledge, if any, of the periodicity for the par- 
ticular signal of interest, and require an algorithm 
that will allow us to find it. The results of § 3.3 and 
§ 4 will serve as a guide in developing this. First, 
we subtract the mean from the signal and perform 
the continuous wavelet transform given by Equa- 
tion (1), only multiplying by a factor of y/Si to 
allow proper normalization and replacing the in- 
tegral with a sum. Then, we divide the continuous 
coefficients by the standard deviation of the time 
series to assure that the expectation value of |/| 2 
for white noise is unity and cross them with the 
coefficients for a number of mock (analyzing) sig- 
nals, each of which has a period r„ = r„_i + St. 
For each of these, we average over the relevant 
region in wavelet space as in Equation (29), and 
look for extrema. The extrema of F c (t„) then cor- 
respond to characteristic time scales, related by 
r a = 0.973f„, where f n is the period of the mock 
signal corresponding to the extremum. 

The nth mock signal is given by 

/„ (t) = A n cos (yt + <f^j , (34) 

where we use the subscript n instead of m on the 
parameters to emphasize that they are for the n th 
mock signal. We choose a value of A n such that 
the global wavelet spectrum of the analyzing sig- 
nal is normalized so that its maximum value is 
unity This provides a convenient normalization 
among the analyzing signals, assuring that they 
all have the same power. In this analysis, 100 
mock signals are used to analyze the actual time 
series, resulting in a step size of St ~ 0.045 for 
most sources (this varies slightly because not all 
of the sources have the same time window). We 
set 4> n equal to zero, as |/ c | 2 is independent of the 
phase difference between the signals. Many of the 
signals analyzed in this paper have typical values 
for the time range on the order of a — 80 yr and 
b = 101 yr, which arise from using the year 1900 
as a baseline. The dilation values are chosen such 



that 2St < I < NdataSt/3 yr, which allows us to 
be conservative when omitting edge effects and to 
admit the approximations used in Equations (8) 
and (13). 

When plotting the coefficients, those that fell 
within the cone of influence at any given dilation 
were set to a constant value, thus masking the 
edge effects. Values of lj are shown logarithmi- 
cally so as to allow more sensitivity in the dilation 
region corresponding to shorter mock periods. It 
is helpful to examine the global wavelet spectrum 
of the time series (Eq. [20] ) , and to compare with 
the cross-wavelet results, so we include it in our 
analysis. We assume a background spectrum Pj of 
the form AR(1) (Eq.[18]), and estimate the lag-1 
autocorrelation coefficient by calculating the lag-1 
and lag-2 autocorrelations, cti and a 2 . The lag- 
1 autocorrelation coefficient is then estimated as 
a = (cti + ^/o2)/2. The background spectrum Pj 
then allows us to compute confidence levels using 
the results of § 4. We plot contours on the continu- 
ous and cross- wavelet transforms corresponding to 
the 90% confidence level for an AR(1) process. We 
also plot 90% and 75% correlated noise confidence 
levels on the plots of the global wavelet spectrum 
f\ G and the averaged cross- wavelet power F c (r„), 
as well as 90% confidence levels for white noise 
(Pj = 1). It should be noted that by confidence 
level, we mean that this given percentage of as- 
sumed background (i.e., noise) processes will pro- 
duce behavior less than what is seen in the wavelet 
transforms. For the time series analyzed here, the 
90% confidence levels for correlated noise mark 
where 90% of AR(1) processes will produce values 
less than the 90% values, and likewise for the 75% 
confidence levels and the 90% confidence levels for 
white noise. 

To satisfy the admissibility condition (4) we 
choose a value of k — 6. In reality, Equation (4) 
is not formally satisfied for the Morlet Wavelet of 
Equation (1). However if k > 5, the admissibil- 
ity condition is satisfied to within the accuracy of 
computer algorithms using single precision arith- 
metic. 

5.2. Simulations 

To test our technique, we ran several simula- 
tions. First, to investigate the ability of our tech- 
nique to detect a purely periodic signal in the pres- 
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ence of noise we used a signal of the form 

/„(*) = 2 sin(7rf + 1.1) +JVn(t), (35) 

where TV is an amplitude factor and n(t) a noise 
function that generates white noise. We show the 
continuous and cross transforms for N — 10 in Fig- 
ure 3. The cross wavelet plot is that for a mock 
signal of period t„ = 2.059 yr. Plots of the global 
wavelet spectrum and F c for this simulation are 
shown in Figure 5a. As can be seen in Figure 5a, 
the period is easily recovered using the technique 
with noise of amplitude five-times that of the am- 
plitude of the sinusoid (signal-to-noise of 0.2). Be- 
cause signal-to- noise values as low as S/N = 0.2 
are never realized, our technique is is valid for the 
time series studied here. 

In addition, we also applied the algorithm to 
a sinusoid that changes amplitude and phase over 
the time window. The transforms are shown in 
Figure 4, and the power spectra in Figure 5b. A 
change in amplitude affects the amplitude of the 
wavelet coefficients, as would be expected. The ef- 
fects of a discontinuous change in the phase <j) a can 
be seen as disturbing the structure of the trans- 
forms in these regions. The technique gives a char- 
acteristic period of 2.0 yr, showing that a sudden 
change in the phase of a signal does not effect the 
technique. 

In practice, we are using this technique to inves- 
tigate quasiperiodic behavior and to give a char- 
acteristic time scale for a time series, rather than 
attempting to recover a periodic signal buried be- 
neath noise (although this certainly may be done). 
To illustrate our technique in the case of pure 
noise, we applied the cross- wavelet technique to 
Gaussian white noise and correlated noise. The 
wavelet transforms for a white noise signal are 
shown in Figure 6. As can be seen, temporary 
quasiperiodic structure can arise for white noise, 
however the signal is clearly distinguished from 
more coherent signals by the complex structure 
throughout the plot, particularly in the low dila- 
tion (high frequency) region. Plots of the global 
wavelet power spectrum and the averaged cross 
power F c are shown in Figure 8a. Although there 
are distinct peaks in the global spectrum at small 
scales, the average power of the cross- wavelet, F c , 
does not show distinct peaks at these scales. We 
interpret this as being the result of smoothing in 
scale as well as time when calculating F c ; the 
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Fig. 2. — Global wavelet power spectrum f\ G {lj) 
and the average cross power F c (t„) for a sinusoid 
of period r a = 2.0 yr. Also shown are the confi- 
dence levels for white and correlated noise, assum- 
ing a lag-1 autocorrelation of a = 0.911362. 
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Fig. 5. — Global power spectrum f\ G {lj) and av- 
erage cross power F c (t„) for the sinusoids of Fig- 
ure 3 (top two panels) and Figure 4 (bottom two 
panels). Confidence levels for correlated and white 
noise arc shown. 
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global wavelet spectrum is averaged in time for 
a given dilation lj, and can result in peaks even 
in the case of non-stationary power, as is seen in 
the continuous plot for white noise. However, by 
crossing with an analyzing signal, and averaging 
in dilation, we can smooth out spurious extrema, 
as is seen in the plot of F c . 

To illustrate the technique for correlated noise, 
we used a test signal of the type AR(1) with lag-1 
autocorrelation a = 0.9. Many of the PR sources 
had lag-1 autocorrelations of a <~ 0.9, and we 
found it helpful to simulate an AR(1) process of 
a similar to what we infer for the data. Figure 7 
shows the transforms for this signal. As expected, 
there is more activity at larger scales (lower fre- 
quencies). Coherent activity does appear in the 
case of this type of correlated noise, however the 
power is not stationary and can not be character- 
ized by a particular scale. The global power spec- 
trum and F c can be seen in Figure 8b. Peaks ap- 
pear at similar periods for both, however the plot 
of F c shows that the peaks are less significant than 
those seen in the global power spectrum. The rea- 
son for this can be seen in the continuous plot; the 
power appears to split into two time scales around 
t' = 90. When we cross the continuous transform 
with our analyzing signal transform, and average 
over both translation and dilation, we see that the 
time scale of ~ 3.8 yr is not as significant as it 
appears in the global wavelet spectrum. 

We can see from the simulations, as well as 
from the results of applying the transforms to the 
PR sources, that the average cross power, F c {r n ), 
places more emphasis on the smaller scales and pe- 
riods than the global power spectrum. This is to 
be expected, as the amplitude of the global power 
spectrum for the time series of interest, f\c(lj), 
has a dependence on the scale lj (see Eq.[8]). If 
we do not normalize the global power spectrum 
of the mock signal, |/ m | 2 7 such that its maximum 
value is unity for all periods r„, then the ampli- 
tude of | / m | 2 will also have a scale dependence. 
This would result in the amplitude of F c being 
dependent on the scale, lj as well (note that F c 
would not be dependent on lj because we divide 
by the scale in Eq.[29]). However, normalizing the 
mock power spectrum |/ m | 2 so that it its maxi- 
mum value is unity removes the scale dependence 
in |/ m | 2 , and thus in F c . 

In addition, we applied the cross-wavelet trans- 
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Fig. 8. — Global power spectrum f\a{lj) and av- 
erage cross power F c (r n ) for the white noise signal 
of Figure 6 (top two panels) and the correlated 
noise signal of Figure 7 (bottom two panels). Con- 
fidence levels for correlated and white noise are 
shown. 
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form technique to the BL Lacs OJ 287 and AO 
0235+164, as previous analysis has shown evi- 
dence for periodicity in these sources. For OJ 287 
the technique detected a characteristic time scale 
of <~ 1.6 yr, with a shorter time scale of <~ 1.1 yr, 
confirming results found earlier by Hughes, Allcr, 
& Aller (1998). The cross- wavelet analysis also 
confirmed earlier results for AO 0235+164, giving 
time scales of <~ 1.9 and <~ 3.3 yr, which is in 
good agreement with the 3.61 yr found by Roy et 
al. (2000) using a Lomb-Scargle periodogram. We 
were not able to confirm the longer periods found 
by Roy ct al. (2000), as they fell outside of our con- 
dition of at least four cycles in the time-window. 

6. THE DATA ANALYSIS 

The results of using the algorithm described 
in § 5.1 to search UMRAO data on the Pearson- 
Readhead survey sources are given in Table 1. Ob- 
servations that we considered to have insufficient 
data to give a reliable analysis of quasiperiodicity 
(about half) are denoted by D. Such observations 
were excluded because we required at least 100 
data points for a given time series, as we assumed 
that 100 points are needed to adequately define 
its character. When analyzing the light curves, 
we used the entire time window, which varied be- 
tween observing frequencies and sources. Since we 
are using a finite step size St <~ 0.045 yr, we do 
not expect to use an analyzing signal with period 
T n = T a , but rather with period t„ s=a r a . However, 
because we are analyzing quasiperiodic behavior 
we do not expect the sources to have a well-defined 
period (i.e., a spike in Fourier space), but rather a 
characteristic time scale. The error resulting from 
the finite step size of 0.045 yr is irrelevant in this 
context, as there is no sense found in assigning an 
extremely precise value for a quasiperiodic source. 
We only record time scales at or above the 90% 
confidence level. 

Many of the sources that exhibited quasiperi- 
odic variations had characteristic time scales be- 
tween one and four years, with an average time 
scale for all quasiperiodic sources of 2.4 yr and a 
standard deviation of 1.3 yr. This, of course, is 
to be expected as our observing interval is around 
twenty years, and we are only interested in varia- 
tions with at least four possible repetitions across 
the interval. Table 2 shows statistics with respect 



to source type and Figure 9 shows histograms for 
the number of characteristic periods found with re- 
spect to source type, binned every 0.25 yr. Out of 
thirty total sources with sufficient data, eighteen 
were quasars, four were galaxies, and eight were 
BL Lacs. A little over half of the sources showed 
evidence for quasiperiodicity in at least one ob- 
serving frequency. Column six shows the ratio 
of recorded time scales to individual time series, 
which gives an idea of the average number of time 
scales per source. For instance, if three quasars 
had sufficient data, and these quasars had eight 
individual time series between them over all the 
observing frequencies, with five of these individ- 
ual time series exhibiting QPVs, then the ratio of 
recorded characteristic periods to individual time 
series would be 5/8 for quasars. A value of 0.33333 
would mean that on average most of the sources 
for this particular type showed evidence for one 
characteristic time scale over the three observing 
frequencies, and two time scales over three observ- 
ing frequencies for a value of 0.66667. As can be 
seen from the table, on average there was about 
a little more than one characteristic period ob- 
served for every three observing frequencies for all 
sources. The last column gives the ratio of sources 
with quasiperiodic behavior in more than one ob- 
serving frequency to the total number of quasiperi- 
odic sources. Roughly half of the quasiperiodic 
sources exhibit quasiperiodicity across more than 
one observing interval. It appears that quasars 
exhibit characteristic time scales more frequently 
than BL Lacs, as about two-thirds of the QSOs 
showed QPVs, while only about on-third of the 
BL Lacs did; however, all of the BL Lacs showing 
QPVs repeated the behavior in more than one ob- 
serving frequency, while only half of the quasars 
did. This may well reflect the fact that the BL 
Lacs of the PR survey sources generally have flat 
spectra (Allcr ct al. 2002). In addition, there 
does not appear to be a distinction among opti- 
cal classes regarding the length of the time scales. 
Among sources showing QPVs in more than one 
observing frequency, the characteristic time scales 
were similar among the different observing fre- 
quencies. However, it is difficult to perform a com- 
plete statistical analysis of quasiperiodicity with 
respect to source type, as only 30 sources were 
analyzed. In addition, there does not appear to 
be a distinct relationship between time scales and 
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observing frequency; characteristic time scales do 
not consistently lengthen with increasing observ- 
ing frequency or vice versa. The plots of one par- 
ticularly promising source for periodic behavior, 
the quasar 0804+499, can be seen in Figures 10 
and 11 at an observing frequency of 4.8 GHz. 

Comparing the results for the PR sources with 
the results expected for correlated noise (§ 4), we 
find it unlikely that most of the PR sources with 
characteristic time scales are exhibiting behavior 
that can be attributed to an AR(1) process. If 
the time series for the PR sources were the result 
of AR(1) processes, we would expect the ratio of 
recorded time scales to time series, Qt/T, to be 
~ 0.1, as we record characteristic time scales at or 
above the 90% confidence level. However, as can 
be seen in Table 2, Q t /T = 0.42. 

The lack of uniform quasiperiodic behavior 
across observing frequencies is unexpected for 
some of the sources with flat spectra. We have 
explored the light curves for these sources and 
conclude that, while generally flat, they do exhibit 
behavior that is different across the three frequen- 
cies. Localized activity is sufficient to remove or 
significantly weaken evidence for quasiperiodic be- 
havior during that time interval, and if quasiperi- 
odic behavior was noted at earlier or later times, 
it often did not span enough cycles (around four) 
to warrant recording a characteristic time scale. 
Typically such events are not sufficient to affect 
the average spectral index, however they can result 
in changes in the small-scale structure of the light 
curves which is detected by the wavelet technique. 
As an example, Figure 12 shows the continuous 
transforms for 1823+568 at 8.0 GHz and 14.5 
GHz, and the wavelet power spectra are shown 
in Figure 13. Although the light curves appear 
similar, the continuous transforms show behavior 
that is different between the two observing fre- 
quencies. Comparison with the power spectra and 
the continuous plots explains the discrepancy; the 
transform for the 8.0 GHz curve exhibits a tran- 
sient time scale of ~ 2.0 yr, seen after <~ 1993, 
whereas the 14.5 GHz curve exhibits a longer time 
scale of ~ 3.3 yr that persists throughout the time 
window. We address the lack of common char- 
acteristic time scales with respect to observing 
frequency further in § 7. 

In addition, we also note that for the 14.5 GHz 
data for 1823+568 the global wavelet spectrum 
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Fig. 9. — Histograms for the quasipcriodicity 
statistics. Time scales are binned every 0.25 years. 
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Fig. 11. — Global power spectrum f\c{h) an d 
the average cross power F c (r n ) for 0804+499 at 
4.8 GHz. The three peaks corresponding to the 
three time scales seen in Figure 10a are seen in 
both plots. The <~ 2.7 yr time scale is the most 
prevalent, as expected from visual inspection of 
the continuous plot. Also shown are the confidence 
levels for correlated and white noise. 
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Fig. 13. — Global power spectrum f\a(lj) and 
the average cross power F c (t„) for 1823+568 at 
8.0 GHz (top two panels) and 14.5 GHz (bottom 
two panels). Although the <~ 2.0 yr time scale ac- 
tivity seen post-1990 in the continuous transform 
for the 8.0 GHz data is significant (see Figure 12a), 
these plots show that the source docs not globally 
show any behavior in 8.0 GHz that is not incon- 
sistent with an AR(1) process. For the 14.5 GHz 
data, power is shared over a broad range in scale in 
both f\ G and F c (t„), however the cross transform 
smooths the two peaks seen in the global spectrum 
into one distinct peak. The single time scale inter- 
pretation given by the cross technique is visually 
reinforced in the continuous plot (see Figure 12b). 



f\c(h) an d the average cross power F c (r n ) differ 
in assigning a time scale. The global spectrum re- 
veals evidence for two time scales, one at ~ 2.9 yr 
and the other at ~ 3.8 yr, while the average cross 
power shows evidence for a time scale at <~ 3.3 
yr. Comparison between the two leads us to con- 
clude that the characteristic time scale seen in the 
plot of F c is more accurate. The global wavelet 
spectrum shows a broad spread of power between 
the two time scales, where the time scales are in- 
ferred from the somewhat poorly resolved peaks. 
The average cross power plot also shows a broad 
spread of power, but with only one peak. This be- 
havior is similar to what was described in § 5.2 for 
the case of white noise; because the cross-wavelet 
technique smooths in dilation as well as in trans- 
lation, and acquires information regarding time 
scales from the analyzing time series, it smooths 
the two poorly resolved peaks in the global wavelet 
spectrum into one broad peak. Analyzing the con- 
tinuous plot, we see that the plot of F c better rep- 
resents the behavior seen in the continuous trans- 
form, as power appears to shared among a range 
of scales of ~ 1.0 yr, rather than two distinct time 
scales separated by about a year. It is more ap- 
propriate to assign one characteristic time scale 
in this case rather than two. Behavior similar to 
what is seen in the 14.5 GHz data was seen in the 
4.8 GHz data as well, but with less significance. 

A comparison of the results from Table 1, with 
those from the structure function analysis per- 
formed by Hughes, Allcr, & Aller (1992) finds 
broad agreement, in the sense that most of the 
sources that we find to have characteristic time 
scales also exhibited time scales in the structure 
function analysis, and likewise for those sources 
lacking time scales. After comparing the time 
scales from the structure function analysis for 
those sources that had time scales short enough to 
meet our requirement of four repetitions over the 
time window (i.e., shorter than 4-5 years for most 
sources), we find that many of the time scales we 
find are comparable to those found from the struc- 
ture function analysis. For instance, we deduce a 
time scale of <~ 2.6 yr for the source 0836+710 at 
4.8 GHz, while the structure function finds that 
the time series for this source is not correlated 
above a time scale of 2.88 yr; also, we find a 
characteristic time scale of ~ 2.0 yr for the BL 
Lac 1803+784, and the structure function analysis 
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gives a time scale of 1.86 yr. Although the struc- 
ture function analysis gives a measure of the time 
scale above which variations appear to be uncorre- 
cted, which is not the same as the time scale that 
we measure, we find that the results found from 
the cross- wavelet technique agree with those found 
the earlier structure function analysis of Hughes, 
Aller, & Aller (1992). 

Color postscript plots for the continuous and 
cross- wavelet transforms, as well as the plots for 
the global wavelet spectrum and the average cross 
power, are available for all sources from the UM- 
RAO website at 

http://www.astro.lsa.umich.edu/ 
obs / radiotel / prcwdata. html . 

7. CONCLUSIONS 

We conclude that, complimented with the con- 
tinuous wavelet transform, the cross-wavelet tech- 
nique can be an effective tool in the search for 
quasipcriodicity of a time series. Comparing the 
results for the Pearson-Readhead survey sources 
with that expected for a correlated noise process 
of the form AR(1), we conclude that the observed 
quasiperiodic behavior is unlikely to be the result 
of an AR(1) process, as about 40% of the time 
series for the PR sources had characteristic time 
scales compared to about 10% expected for an 
AR(1) process. The observed time scales may be 
the result of a type of correlated noise that is not 
AR(1), and it certainly may not even be station- 
ary noise. However, even if the QPVs arise from 
correlated noise, it is meaningful to explore the 
characteristic time scale of a time series, and such 
results provide a useful diagnostic of the underly- 
ing variations. 

After applying the cross-wavelet algorithm to 
the Pearson-Readhead VLBI survey sources, anal- 
ysis revealed evidence for quasi-periodic variations 
in <~ 57% of the sources, as well as evidence that 
- 67% of quasars, - 38% of the BL Lacs, and 50% 
of the galaxies have quasiperiodic behavior in at 
least one observing frequency. The sources were 
observed to have a mean characteristic time scale 
of 2.4 yr, with standard deviation of 1.3 yr. 

Because the analyzed radio band variations for 
the PR sources originate in parsec-scale jets, per- 
turbations that propagate at c will lead to observ- 
able fluctuations if they span the flow, guarantee- 



ing time scales of order years. In addition, coher- 
ent perturbations that arise from the excitation 
of certain modes of oscillation of the flow, which 
could give rise to quasiperiodic behavior with time 
scales comparable to the dynamical response time 
of the flow (i.e., years), have been shown to arise 
quite naturally (Hardee et al. 2001). We conclude 
that our results are in good agreement with the 
characteristic time scales that we would expect to 
observe based on the nature of these objects. 

There is a transition region where the jet 
changes from optically thick to optically thin, with 
optical depth r = 1 (Cawthorne 1991). Because 
this region varies with the observing frequency, we 
are looking at a different physical location in the 
jet at each observing frequency. It is likely that 
the time scale of quasiperiodic variations is de- 
pendent on their location in the jet. Naively, we 
would expect to probe larger scale regions with 
longer characteristic time scales, the lower the 
observing frequency chosen. The results of our 
analysis certainly find no such correlation of time 
scale and frequency. However, jets may accel- 
erate due to adiabatic expansion, or decelerate 
due to entrainment, leading to a change in bulk 
Lorentz factor with position. In addition, cur- 
vature is now known to be a common feature of 
these flows. A change in flow speed and/or flow 
orientation with respect to the observer can lead 
to a significant change in Doppler factor, and thus 
to the observed time scale of activity, masking any 
simple trends. Furthermore, local jet properties 
and ambient conditions play a major role in de- 
termining what modes of the flow exist, and may 
be excited, and observations that probe different 
physical scales might well reveal activity in one 
frequency band but not in another. 

Although the number of sources analyzed here 
is in no way exhaustive, we see no reason why 
such quasiperiodic behavior should be confined to 
the Pearson-Readhead survey sources, and we find 
it likely that many active galactic nuclei exhibit 
quasiperiodicity. Only four galaxies had sufficient 
data for analysis, and so are tabulated but not 
discussed. A statistical analysis on the results of 
applying the cross-wavelet technique to a greater 
number of sources would lead to a more interest- 
ing and conclusive comparison of quasiperiodic be- 
havior and source type, as well as evidence for 
quasiperiodicity in active galactic nuclei in gen- 
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eral. In addition, the cross-wavelet transform may Torrence, C, & Compo, G.P. 1998, Bull. Amer. 

also be a useful tool for analyzing the correlation of Meteor. Soc, 79, 61 

a source between different observing frequencies. 

We would like to thank the anonymous ref- 
eree for many helpful and insightful comments 
that have contributed to significant improvement 
of this manuscript. 
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Table 1 

Characteristic periods found using the cross-wavelet algorithm on the 



Pearson-Readhead VLBI survey sources. All periods are in 
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BL LAC 


BL 


1.4 


856 


3.7 


1441 


3.5, 1.6, 0.7 


1243 


2229+391 


G 


D 




D 




D 




3C 452 


G 


D 




D 




D 




2351+456 


Q 




148 


D 




3.8 


167 


DA 611 


G 


D 




D 




D 





Note. — D signifies insufficient data, - signifies no detected quasipcriodicity. The number of data 
points are given only for those sources with sufficient data. If two time scales are given, the more 
prominent one is listed first. 
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Fig. 1. — (a) Continuous wavelet transform for a 
sinusoid of period r a = 2.0 yr. (b) Cross-wavelet 
transform for the sinusoid with another sinusoid 
of period r m = 2.059 yr. The contours on both 
power (modulus) plots are for the 90% confidence 
level for correlated noise, assuming a lag-1 auto- 
correlation of a = 0.911362. 



Fig. 3. — (a) Continuous wavelet transform for a 
sinusoid of period r a = 2.0 yr with white noise 
added, S/N = 0.2. (b) Cross-wavelet for the sinu- 
soid of (a). The contours on both power plots are 
for the 90% confidence level for correlated noise, 
assuming a lag-1 autocorrelation of a = 0.283620. 



Fig. 4. — (a) Continuous wavelet transform for a 
sinusoid of period r a = 2.0 yr and varying ampli- 
tude and phase, (b) Cross-wavelet transform for 
the sinusoid of (a). The contours on both power 
plots are for the 90% confidence level for corre- 
lated noise, assuming a lag-1 autocorrelation of 
a = 0.894581. 



Fig. 10. — (a) Continuous wavelet transform for 
0804+499 at 4.8 GHz. The periodic behavior is 
easily seen as repetition in the real plot. The mod- 
ulus plot shows the significance of the behavior; 
particularly notable is the ~ 2.7 yr period, (b) 
Cross-wavelet transform for 0804+499 at 4.8 GHz 
with an analyzing signal of period r n = 2.86 yr. 
The data is well correlated with a sinusoid of this 
period, as is evidenced by the uniformity in the 
cross-wavelet plots. The contours on both power 
plots are for the 90% confidence level for corre- 
lated noise, assuming a lag-1 autocorrelation of 
a = 0.913. 



Fig. 6. — (a) Continuous wavelet transform for a 
Gaussian white noise signal. White noise is eas- 
ily recognizable by its erratic and complex struc- 
ture in wavelet space, particularly at small dila- 
tions, (b) Cross-wavelet transform for the Gaus- 
sian white noise signal, corresponding to an ana- 
lyzing signal of period r„ = 0.736 yr. The contours 
on both power plots are for the 90% confidence 
level for white noise. 



Fig. 7. — (a) Continuous wavelet transform for a 
correlated noise signal of the form AR(1), with lag- 
1 autocorrelation coefficient a = 0.9. (b) Cross- 
wavelet transform for the correlated noise signal. 
The signal is crossed with a mock signal of period 
t„ = 3.901 yr. The contours on both power plots 
are for the 90% confidence level for an AR(1) pro- 
cess of a — 0.9. 



Fig. 12. — (a) Continuous wavelet transform for 
the BL Lac 1823+568 at 8.0 GHz. The contours 
are the 90% confidence level for an AR(1) process 
of a = 0.864. Note the significant behavior seen 
after 1990. (b) Continuous wavelet transform for 
1823+568 at 14.5 GHz. The contours are the 90% 
confidence level for an AR(1) process of a = 0.987. 
One can see that power is distributed over a broad 
range in scale for this plot, centered about a dila- 
tion corresponding to a time scale of ~ 3.3 yr. Al- 
though the two time series share similarities, there 
are enough differences to account for the varying 
behavior seen in the transforms. 
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Table 2 
quasiperiodicity statistics. 



Type 


Sources 


Time Scries 


Mean Time Scale (yr) 


Stnd. Dev. (yr) 


Q/s 


Qt/T 


M/Q 


1 


2 


3 


4 


5 


6 


7 


8 


Quasars 


18 


51 


2.3 


1.4 


0.67 


0.43 


0.5 


Galaxies 


4 


10 


2.0 


1.1 


0.5 


0.4 


0.5 


BL Lacs 


8 


24 


2.3 


1.1 


0.38 


0.42 


1.0 


All Types 


30 


85 


2.4 


1.3 


0.57 


0.42 


0.59 



Note. — Q/S is the ratio of sources exhibiting quasiperiodicity in at least one observing frequency to the total 
number of sources. Qt/T is the ratio of the total number of time scales to the total number of time series. M/Q 
is the ratio of sources that showed quasiperiodicity in more than one observing frequency to the total number of 
quasiperiodic sources. 
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